#if 0
            // midpoint rule; psi_na=psi_nb=0.5 at midpoint of segment
            sint = 0.5 * user->gN_fcn(0.5*(ax[na]+ax[nb]),0.5*(ay[na]+ay[nb])) * ls;
            if (abfn[na] != 2)  // node at end of segment could be Dirichlet
                aF[na] -= sint;
            if (abfn[nb] != 2)
                aF[nb] -= sint;
#else
            // two-point gauss rule
            const double z0 = (1.0 - 1.0/sqrt(3.0)) / 2.0,
                         z1 = (1.0 + 1.0/sqrt(3.0)) / 2.0;
            double       x0, x1, y0, y1, gN0, gN1;
            x0 = z1 * ax[na] + z0 * ax[nb];
            x1 = z0 * ax[na] + z1 * ax[nb];
            y0 = z1 * ay[na] + z0 * ay[nb];
            y1 = z0 * ay[na] + z1 * ay[nb];
            gN0 = user->gN_fcn(x0,y0);
            gN1 = user->gN_fcn(x1,y1);
            if (abfn[na] != 2) { // node at end of segment could be Dirichlet
                sint = 0.5 * (gN0 * z1 + gN1 * z0) * ls;
                aF[na] -= sint;
            }
            if (abfn[nb] != 2) {
                sint = 0.5 * (gN0 * z0 + gN1 * z1) * ls;
                aF[nb] -= sint;
            }
#endif
